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Abstract 

A kinetic description of lattice-gas automaton models for reaction-diffusion 
systems is presented. It provides corrections to the mean-field rate equations 
in the diffusion-limited regime. When applied to the two-species Maginu 
model, the theory gives an excellent quantitative prediction of the effect of 
slow diffusion on the periodic oscillations of the average concentrations in a 
spatially homogeneous state. 
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I. INTRODUCTION 



In this paper we describe how a relatively simple theory quantitatively explains the 
deviations from mean-field behavior that occur in diffusion-limited chemical reactions. The 
modeling of chemical reactions in spatially extended systems is an interesting application 
of a class of microscopic models called "lattice-gas automata" [||. Space, velocity, and time 
are all discrete in such models, which simplifies implementation on computers as well as 
theoretical analysis. Lattice-gas automata (LGA) provide a fiexible tool for studying the 
various phenomena resulting from the interplay between reaction and diffusion 0. 

Here we will not be concerned with chemical pattern formation, but instead we will 
consider a spatially extended two-species model exhibiting coupled periodic oscillations of 
the concentrations of both species in a spatially homogeneous state. If the reactions are 
slow compared to the diffusion, then mean-field or Boltzmann theory equations give an 
excellent description of the reaction kinetics. This is the so-called reaction-limited regime. 
In the opposite diffusion-limited case however, when the diffusion is slow compared to the 
reactions, there is no time to equilibrate after a reaction before another reaction occurs. 
Consequently, equal-time correlations will be present that invalidate the Stosszahlansatz 
or molecular chaos assumption used to derive the Boltzmann equation. Therefore in the 
diffusion-limited regime the behavior of the system is seriously modified. 

A condition that guarantees the absence of correlations in the equilibrium state of lattice- 
gas automata is the so-called detailed balance (DB) condition. Reactive LGA's in the 
diffusion-limited regime violate DB. A systematic theory for LGA's violating DB has recently 
been developed by Ernst and coworkers P,^. In the present paper we apply this theory to 
calculate corrections to the Boltzmann equation. A similar method has been developed by 
Boghosian and Taylor PJ^. 

The organization of this paper is as follows. In section |I| we define the model used. We 
present the ring kinetic theory in section |ITT| , and compare it with computer simulations in 
section We end with a discussion in section |V[ 



II. THE MODEL 

A. Reactive Lattice-Gas Automaton 

In a lattice gas automaton particles live on a regular lattice, C, so that their positions can 
only take a a limited set of values corresponding to the nodes of the lattice. The velocities 
are also restricted, and must be equal to unit vectors oriented along the the links connecting 
the neighboring nodes. We denote this set by {cj; I < i < b} where b is the coordination 
number of the lattice. The square lattice, with 6 = 4, will be used in this paper as it has 
sufficient symmetry to properly describe the diffusive problem that we are considering. We 
further impose an exclusion principle requiring that no more that one particle can be at 
the same node with the same velocity. As a consequence there can be at most b particles 
per node, i.e., one per link. The state of the LGA is fully described by a set of boolean 
occupation numbers {si(r); I < i < b, r G £}, where Si{r) equals 1 if there is a particle at 
node r with velocity Cj and otherwise. 
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For multi-species models with m types of reactants, such as the Maginu model where 
m = 2, we have to introduce different types of particles. The exclusion principle has to 
be modified in order to allow for the coexistence of several species. We adopt the coupled- 
lattice model described in 0. In this approach particles of different types live on separate 
lattices, and only interact when a chemical reaction occurs. The exclusion principle is applied 
independently to each lattice. However, for the sake of compactness in the mathematical 
derivations, we can extend the former set of occupation numbers Sj(r) to a new set {sj(r); 1 < 
i < mb, r G £}, in such a way that channels 1 < i < b are reserved for particles of species 
1, channels b + 1 < i <2b for species 2, etc. The number of particles of type p is given by 



pb 



ap[r) 



E ^^W- (1) 



i=l+(p-l)6 



A time evolution step is the composition of two substeps, defined as follows. First, at 
each node independently a reactive collision takes place, during which a pre-reaction state 
s(r) = {sj(r), I < i < mb} is replaced by a post-reaction state cr(r) in a stochastic process 
governed by a set of transition probabilities Ag^. The reactive collision is followed by a 
propagation step, during which all particles are moved to neighboring nodes r -|- Cj in the 
direction of their velocities. 

Let us describe the reactive collision step in detail. The chemical reaction we want to 
simulate is described by 

aiXi + 02X2 + ... + amXm — > (3iXi + (32X2 + . . . + (2) 

and occurs at a rate P{a,(3), where cx = {ai,a2, ■ ■ ■) and (3 = {Pi, P2, ■ ■ ■) specify the 
number of particles before and after reaction, and Xp represents species p. The outcome of 
the chemical reaction only depends on the number of particles of each species, {ap{s); 1 < 
P < fn}, present at the node before the reaction, not on the velocity distribution. After 
the reaction, the jSp particles of each species are randomly redistributed over the b available 
velocity directions (this random redistribution models the diffusion process), which can be 
done in 5!/(/3p)!(6 — /?p)! ways for species p. Thus, the transition probability from precoUision 
state s to postcoUision state a is given by 



A.. 



- {(3p{a))\{b - Pp{a))\ 



P{cx,f3). (3) 



Note that the normalization I]^ ^scr = 1 follows from the normalization J^p P{f^, ^) = 1- 



B. Maginu Model 

The Maginu model is a two species model that exhibits a variety of behavior. It is 
described by the following equations for the concentrations x and y of the two species : 

_ = a; - x^S -y + D,.V\ 

^ = (x-ky)/c + DyV'y (4) 
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with c > and < k < 1. The constants and Dy are the diffusion coefficients for 
the two species respectively. Depending on the parameters, the model can exhibit Turing 
structures (when is very different from Dy) as well as periodic behavior. Here we will 
solely be interested in the case = Dy, where the system develops a stable limit cycle in 
a homogeneous state. This limit cycle shrinks as the chemical reaction rate increases. 

The Maginu model as defined by Eq. (^) is not directly useful since the concentrations 
can become negative, and therefore cannot be simulated with an LGA |^. This problem is 
however easily solved by using the linear transformation 

x=l + x/^12{l + k)/k, 

y = I + yk/ ^12(1 + k)/k (5) 

where x and y are the concentrations of the two species X and Y that we will study. 

Next we have to determine a set of reaction rates P{cx.,l3) for the LGA that gives rise 
to the macroscopic behavior defined by Eqs. (|^) and (|^). The matrix P{a,f3) is needed in 
numerical simulations as well as in the theory presented in the next section. In Ref. a 
method for constructing P{cx.,f3) has been extensively discussed, and we will not give the 
details here. We will however, adopt the rules of Ref. |^, where the number of particles is 
only allowed to change by ±1 during the reaction. The matrix P{a, (3) is then uniquely 
specified. 

An important point in the definition of the collision rules is the introduction of a time 
scaling parameter, s, which allows us to control whether the system is in the reaction- 
limited or in the diffusion-limited regime (see Ref. for details). For large values of s we 
have P{a,f3) ~ 6{a,f3) (where 6 is the Kronecker delta): chemical reactions occur at a 
very slow rate. This is the diffusion-limited regime, where diffusion is able to maintain the 
homogeneity in the system, and where Eqs. @ and (^ are meaningful, as the conditions 
for their derivation are fulfilled. On the other hand, for small values of s chemical reactions 
occur at a much faster rate, and diffusion is no longer able to maintain spatial homogeneity. 
This is the reaction-limited regime. In the next section we present a theory that explains 
the behavior of the system throughout both regimes. 



III. RING KINETIC THEORY 

In mean-field or Boltzmann approximation all correlations between occupation numbers 
are neglected, and the state of the system is completely specified by the average occupation 
numbers, 

Mr,t) = {s,{r,t)). (6) 
The time evolution of fi{r,t) is given by the nonlinear Boltzmann equation, 

/,(r,t + l) = /,(r,t) + /,[/(r,t)]. (7) 
The nonlinear collision operator is defined as 

Uf] = T.i^i - Si)A,„F{s) = {ai - Si)F. (8) 

S,cr 
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We have introduced {■ ■ ■)f as an average that assumes that the precolhsion state is factorized 
over all channels, so that the probability to find a state s is given by 



i^(^)=n/"(i-/o'"" (9) 

i 

In this approximation, where F{s) is given by Eq. (|^), and the transition rates Asa are those 
of the Maginu model, the nonlinear Boltzmann equation (|^) is equivalent to the mean-field 
rate equations and (^. 

To go beyond the mean-field approximation we consider the pair correlation function, 

C,,(r-r',t) = {6s,ir,t)6s,ir\t)). (10) 

Here we have assumed that the system is translationally invariant. The fiuctuations are 
defined as 6si = Si — fi. A special role is played by the on-node correlations Cjj(0,t); by 
definition the diagonal elements vanish: Cjj(0, t) = 0. We neglect all triplet and higher order 
correlations. In a spatially homogeneous system, where fi{r,t) = fi{t), the time evolution 
of fi{t) is then described by the generalized Boltzmann equation, 

Mt + 1) - Mt) = i,[f{t)] + Y.nM[fit)niit)- (11) 

k<l 

Here the operator /' describes corrections to the Boltzmann collision term /. It is defined 

by 

where gi = ((5sj)^) = /j(l — fi) is the single channel fiuctuation. 

In order to have a complete theory we must provide a time evolution equation for Cjj(r, t). 
To derive this equation we will make the important assumption that the average occupations 
change slowly in time. In fact, as far as the evolution of Cij{r, t) is concerned, we will assume 
that no chemical reactions occur at all, so that the model is purely diffusive. Under this 
assumption, the average occupations in equilibrium are given by 

fP = x (2 = 1,2,3,4) 

fP = y (z = 5,6,7,8), (13) 

where x and y are the average concentrations of species X and Y, respectively. When /j(r, t) 
is close to equilibrium, the approach to equilibrium is given by the linearized Boltzmann 
equation {6fi = fi - /D, 

6f,{r + ci, t + 1) = + n)^,6f,{r, t), (14) 
j 

where Ijj = 6ij and the linearized Boltzmann operator is defined by 

Sfj \ 3i 



n, = ^JM^L,.s.M . (15) 



F 
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Under the assumption of slow reactions we have 
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(16) 



It is natural to assume that fluctuations 6si{r,t) will decay to equilibrium in a manner 
similar to 6fi{r,t), i.e., 



Ssi{r + Ci, t + 1) = ^(U + n)ij6sj{r, t). 



(17) 



However, two fluctuations at the same node will be correlated after collision, even if before 
collision the distribution is completely factorized. This is a consequence of the violation of 
detailed balance 0. The generation of on-node postcollision correlations is quantified by 



(18) 



This expression vanishes in the non- reactive limit s — > oo. The presence of on-node correla- 



tions Cij{0,t) before collision gives rise to corrections to fijj [/(t)], and the full postcollision 
source term is given by (see Ref. [Q) 

i?,,(t) = ^]|°[/(t)] + c,,(o,t) + 5:fi|^,,[/(t)]CH(o,t), 



(19) 



where nf^%[f] = d^nif[f]/dfkdfi. Combining Eqs. and (0) with the definition of 
Cij{r,t) in Eq. (|10D we obtain the ring kinetic equation 



Cij{r + Ci- Cj,t + 1] 



;i - 5r,o) E(l + + ^)jiCki{r,t) + 5r,o5,,[/(t)]. (20) 



kA 



This equation has been derived in a more systematic fashion in Ref. 0. 

The physical interpretation of Eq. ( PU[ ) is as follows. Two fluctuations on the same node 
r that are correlated after collision at time to, will be propagated to neighboring nodes r-|-Cj 
and r + Cj. Due to the collision with other particles at these nodes the correlation will be 
scattered to all directions as described by (1 + ^7). Thus both fluctuations branch into many 
different paths. At time to + ^ the weight of each path is given by the same factor (1/4)"^. 
If two correlated paths end at the same node — a so-called "ring" -collision — they give 
rise to on-node precoUision correlations, C(0,to + t) ~ (l/4)^''i?(to), that change the time 
evolution of the average occupations according to Eq. (|Tl]). The actual value of Cij{r,t) is 
a superposition of "ring" contributions from source terms at all earlier times, although the 
dominant contribution comes from the last few time steps. 

The fact that Eq. (|20D is linear in C allows us to write 
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t-1 

Cij{0, t) = E E ^iMt - t')Bki{t'). (21) 

i'=0 k,l 

Here Kij_ki{t — t') is a memory kernel which does not depend on any of the model parameters 
— although it does depend on the system size L — and thus can be constructed once and 
for all using Eq. (pOl) . This can be done in an efficient manner by exploiting the rotation 
and reflection symmetry of Kij ki{t — t'). 

After an initial fast decay, the memory function decays algebraically for large t, as 
Kij^ki{t) oc t""" with a ~ 1.2 for L = 256. When the ring kinetic theory is evaluated 
numerically, on time scales on the order of 10^ time steps this slow decay leads to the build- 
up of pair correlations that are much larger that what is observed in simulations. This excess 
of correlations would be corrected if we include higher order correlations that are not taken 
into account by the the present form of the ring kinetic theory. Therefore it is desirable to 
cut off the memory kernel for large times, i.e., to set Kij^kiit) = for t > tcutofr- It is natural 
to choose the cut-off equal to the time it takes to travel across the the system: tcutofr = L. 

By rewriting Eq. ( PU] ) in Fourier representation, it can be interpreted in terms of modes at 
different wavevectors q (see Ref. 0). When no reactions occur, the diffusive modes around 
q = and the (spurious) staggered modes around q = (vr, vr) play a special role, since they 
correspond to conserved densities. However, in the presence of reactive collisions there are 
no conserved densities, and all modes are in principle equally important in determining the 
size of the correlations. 

In the next section we will compare theoretical predictions with the results of computer 
simulations. Numerically, the theory of this section is evaluated as follows. At time t = we 
set /i(0) = xo for 1 < i < 4, /i(0) = |/o for 5 < i < 8, and Cjj(r, 0) = 0. To perform a time 
evolution step from time t to time we first use f{t) to calculate the nonlinear Boltzmann 
operator /j[/(t)] and the correction term I'ikiifi^)]- Together with the on-node correlations 
Cjj(0,t) we then use these operators to calculate fiit + 1) with the help of Eq. (|TTp. To 
calculate the evolution of the pair correlation function we use /j(t) and Cjj(0,t) to evaluate 
the source term Bij(t) in Eq. (|I^) and then obtain Cij{r,t + 1) with the help of Eq. (^U|). 
Iteration of the above procedure yields the set {[x(t),y(t)]; t > 0} defining a trajectory in 
the x-y concentration plane. For large times, either a fixed point or a limit cycle is reached. 

IV. COMPARISON WITH SIMULATIONS 

Our simulations were carried out on a 256 x 256 square lattice. The parameters used in 
the simulations were k = 0.9 and c = 2, i.e., identical to those used in Ref. 0. At t = 
the system was prepared in an uncorrelated homogeneous state, with average concentrations 
2^0 = I/O = 0.6. Then we performed the time evolution of the LGA, according to section ||. 
The initial time steps were discarded, as the system needs some time to build up the corre- 
lations that will eventually produce the shrinking of the limit cycle. Once the correlations 
have been created we record the spatially averaged concentration of both species. The scale 
parameter s was varied between s = 2 and s = 20. 

In Fig. ^ the dashed line denotes the limit cycle as it is obtained from the mean-field 
theory defined by Eqs. (^ and (||), assuming that the concentration of both species are 
homogeneous and the term can be neglected. For relatively large values of the time 
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scaling parameter s we expect mean-field theory to be accurate. This is confirmed by the 
simulation data for s = 20, shown as a gray band in Fig. |l|, which are reasonably close to 
the mean-field prediction. The width of the gray band corresponds to the fluctuations in 
the spatially averaged concentrations that occur due to the finite system size. 

When s is decreased correlations become important (measurements show that correla- 
tions are typically 10 times larger for s = 4 than for s = 20) and the diffusion process is 
not able to keep the system homogeneous. As a consequence, different regions in the system 
become desynchronized to a certain degree, and the contribution to the average concentra- 
tion of one region is partially canceled out by out-of-phase contributions from other regions. 
This produces a shrinking of the limit cycle in Fig. ^ as is shown by the simulation data for 
s = 6 and s = 4 in Fig. [l|. The effect is stronger for smaller s. 

It is clear from Fig. |I] that mean-field theory completely fails for the smaller values of s. 
The solid black lines in Fig. |I| represent the limit cycle as it is predicted by the ring kinetic 
theory of section For the values s = 20 and s = 6 shown in Fig. |I| our theory gives 
an excellent quantitative prediction of the shrinking of the limit cycle. For s = 4 there are 
deviations due to higher order effects that are not taken into account. 

Analysis of the ring kinetic theory shows that as s is further decreased, the limit cycle 
shrinks continually, until at s ~ 3 there is an inverse Hopf bifurcation from a limit cycle 
to a fixed point. This bifurcation corresponds to a desynchronization transition, where the 
coherence between different regions is lost completely. It should be noted that this transition 
is of a different character than the Hopf bifurcation that occurs at the mean-field level as a 
function of the model parameters k and c. 

Let us consider the case s = 2 in some detail. Here s is close to the smallest possible 
value (see Ref. 0) and the fluctuations caused by the chemical reactions are strongest here. 
Diffusion is not fast enough to keep the system homogeneous except at very small scales. 
The ring kinetic theory for s = 2 predicts a fixed point located a.t x = y = 1/2. Simulations 
for a system of linear size L = 256 reveal that the average concentrations fluctuate around 
the point x = y = 1/2 in a irregular fashion, and in a range between 0.47 and 0.53. In order 
to assess whether the result of the simulations for s = 2 corresponds to a fixed point or to 
a limit cycle, we compared numerical simulations for three different system sizes: L = 32, 
L = 256, and L = 1024. The concentration of species X versus time is plotted in Fig. |^. 
The vertical scale in all three plots is the same. Clearly, the amplitude of the oscillation 
decreases with the system size. In the L = 32 system the concentration x oscillates with an 
amplitude Sx ~ 0.10; in the L = 256 system we have 6x ~ 0.02, and in the L = 1024 system 
the fluctuations are very small, 6x ~ 0.004. It is reasonable to conclude that for s = 2 the 
correct solution is a stable fixed point, in perfect agreement with our theory. 

In close connection with this last point, we have verified that for s > 4 the limit cycle 
obtained in the simulations is finite and stable, and independent of the size of the system 
up to size L = 1024. For s = 3 our ring kinetic theory predicts a fixed point x = y = 1/2. 
However, simulations are here not conclusive, as systems of intermediate size L = 256 show a 
limit cycle, but large systems do not reach any stationary behavior within available computer 
time. We conclude that the (inverse) Hopf bifurcation from a limit cycle to a fixed point 
at the level of the spatially averaged concentrations in a large enough system must occur 
between s = 2 and s = 4. 

The comparison between ring kinetic theory and simulations has so far been restricted to 
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the shape of the hmit cycle. Figure |l] however does not give any information about the actual 
time evolution of the concentrations, or the period of oscillation around the limit cycle. In 
order to obtain this information, we have plotted in Fig. |^ the average concentration of 
the two species versus time, for both theory and simulations. Figure |^ and B show the 
concentration of particles of type X and F, respectively, for s = 10. Simulation results are 
indicated by a solid line, while the ring kinetic theory is denoted by a dashed line. The 
amplitudes of the oscillation agree quite well, as we already knew from Fig. ^ There is 
however some deviation between theoretical and simulated periods, that causes the curves 
to become slightly out of phase. The difference between both oscillation periods is about 
3%. Figures and D show similar curves for s = 4. Here the agreement is worse, and the 
difference in periods is about 11%. 

Figure^ shows how the oscillation period — normalized by dividing by s — depends on s. 
We have plotted the mean-field value of the period as a dashed line; the ring theory is denoted 
by circles, and simulation results by triangles. It was shown in Fig. |I| that ring kinetic theory 
predicts the shape of the limit cycle quite well down to s ~ 6. It is therefore somewhat 
surprising that the mean- field prediction for the oscillation period, which is s- independent, is 
better than ring kinetic theory for all values of s. To resolve this issue it would be necessary 
to include higher order correlation functions in the theoretical description. This is clearly 
beyond the scope of the present paper. Furthermore, it can be seen that the approach to 
the mean-field period for large s is slow, and even for s = 20 there is a clear deviation of 
about 1%. This effect is probably due to the particular choice of the transition rates, that 
are not able to maintain the local diffusive equilibrium even for high s 0. 

V. DISCUSSION 

In this paper we have shown how a theory that takes into account equal-time pair cor- 
relations provides an excellent explanation of the large deviations from mean-field theory 
observed in diffusion-limited chemical reactions as modeled by lattice-gas automata (LGA). 
Our theory is a straightforward application of the general framework established in the pa- 
pers of Ernst and coworkers [|,|| . It is not restricted to the Maginu model, but is applicable 
to any chemical reaction that can be modeled with an LGA. 

It is in principle possible to include triplet and higher order correlations as well. How- 
ever, the good agreement between theory and simulations indicates that the ring theory of 
section |T| captures the essential physics in a quantitative way. Although the comparison 
between theory and simulations reported here is restricted to the domain of LGA's, we 
expect that mutatis mutandi the general concepts apply equally well to continuous systems. 

We have focused on a particular two-species model exhibiting periodic oscillations of the 
average concentrations. Wu and Kapral |§] have studied a model with more complicated 
temporal behavior — period doubling bifurcations and a transition to a strange attractor, as 
model parameters are changed. They investigated the consequences of spatial fluctuations 
by means of computer simulations. It is an interesting question whether some of the features 
observed in that work can be explained using the theory presented in this paper. As a final 
remark we mention that our theory provides a more microscopic analogue of the Langevin 
equation method used in Ref. [0 to predict the magnitude of spatial density correlations. 
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FIGURES 

FIG. 1. Average concentrations x and y in the Maginu model for several values of the time 
scaling parameter: s = 4, 6, and 20. The outer dashed line corresponds to mean-field theory, given 
by Eqs. (§) and (^). Solid lines correspond to the ring kinetic theory of section |l|. The gray 
bands denote the result of computer simulations performed on square lattices of size 256 x 256; 
their width corresponds to the size of the fluctuations from cycle to cycle. 

FIG. 2. Average concentration x of species X versus time t for the Maginu model at s = 2, and 
system size 32 x 32, 256 x 256, and 1024 x 1024, respectively. The bigger the system, the smaller 
the fluctuations. 

FIG. 3. Concentration of species X and Y versus time t for s = 10 (figures A and B), and for 
s = 4 (figures C and D). Solid lines are the results of the computer simulations (in systems of size 
256 X 256), while the dashed lines correspond to ring kinetic theory. 



FIG. 4. Oscillation period as a function of s. The mean field value is indicated by a dashed 
line. Circles denote the ring kinetic theory prediction of section ED. Triangles are the simulation 
values. Mean-field theory is here in better agreement with simulations than ring kinetic theory. 
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FIGURE 1 
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FIGURE 2 
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